﻿/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms which convert a geocentric set of coordinates to their topocentric equivalent
Created: PJN / 29-12-2003
History: None

Copyright (c) 2003 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/

using System;

namespace AAPlus
{
    public class TopocentricEclipticDetails
    {

        public double Lambda;
        public double Beta;
        public double Semidiameter;
    }

    public static class Parallax
    {
        static double g_AAParallax_C1 = Math.Sin(CoordinateTransformation.DegreesToRadians(CoordinateTransformation.DMSToDegrees(0, 0, 8.794)));


        public static double DistanceToParallax(double Distance)
        {
            double pi = Math.Asin(g_AAParallax_C1 / Distance);
            return CoordinateTransformation.RadiansToDegrees(pi);
        }

        public static double ParallaxToDistance(double Parallax)
        {
            return g_AAParallax_C1 / Math.Sin(CoordinateTransformation.DegreesToRadians(Parallax));
        }

        public static Coordinate2D Equatorial2TopocentricDelta(double Alpha, double Delta, double Distance, double Longitude, double Latitude, double Height, double JD)
        {
            double RhoSinThetaPrime = Globe.RhoSinThetaPrime(Latitude, Height);
            double RhoCosThetaPrime = Globe.RhoCosThetaPrime(Latitude, Height);

            //Calculate the Sidereal time
            double theta = Sidereal.ApparentGreenwichSiderealTime(JD);

            //Convert to radians
            Delta = CoordinateTransformation.DegreesToRadians(Delta);
            double cosDelta = Math.Cos(Delta);

            //Calculate the Parallax
            double pi = Math.Asin(g_AAParallax_C1 / Distance);

            //Calculate the hour angle
            double H = CoordinateTransformation.HoursToRadians(theta - Longitude / 15 - Alpha);
            double cosH = Math.Cos(H);
            double sinH = Math.Sin(H);

            Coordinate2D DeltaTopocentric = new Coordinate2D();
            DeltaTopocentric.X = CoordinateTransformation.RadiansToHours(-pi * RhoCosThetaPrime * sinH / cosDelta);
            DeltaTopocentric.Y = CoordinateTransformation.RadiansToDegrees(-pi * (RhoSinThetaPrime * cosDelta - RhoCosThetaPrime * cosH * Math.Sin(Delta)));
            return DeltaTopocentric;
        }

        public static Coordinate2D Equatorial2Topocentric(double Alpha, double Delta, double Distance, double Longitude, double Latitude, double Height, double JD)
        {
            double RhoSinThetaPrime = Globe.RhoSinThetaPrime(Latitude, Height);
            double RhoCosThetaPrime = Globe.RhoCosThetaPrime(Latitude, Height);

            //Calculate the Sidereal time
            double theta = Sidereal.ApparentGreenwichSiderealTime(JD);

            //Convert to radians
            Delta = CoordinateTransformation.DegreesToRadians(Delta);
            double cosDelta = Math.Cos(Delta);

            //Calculate the Parallax
            double pi = Math.Asin(g_AAParallax_C1 / Distance);
            double sinpi = Math.Sin(pi);

            //Calculate the hour angle
            double H = CoordinateTransformation.HoursToRadians(theta - Longitude / 15 - Alpha);
            double cosH = Math.Cos(H);
            double sinH = Math.Sin(H);

            //Calculate the adjustment in right ascension
            double DeltaAlpha = Math.Atan2(-RhoCosThetaPrime * sinpi * sinH, cosDelta - RhoCosThetaPrime * sinpi * cosH);

            Coordinate2D Topocentric = new Coordinate2D();
            Topocentric.X = CoordinateTransformation.MapTo0To24Range(Alpha + CoordinateTransformation.RadiansToHours(DeltaAlpha));
            Topocentric.Y = CoordinateTransformation.RadiansToDegrees(Math.Atan2((Math.Sin(Delta) - RhoSinThetaPrime * sinpi) * Math.Cos(DeltaAlpha), cosDelta - RhoCosThetaPrime * sinpi * cosH));

            return Topocentric;
        }

        public static TopocentricEclipticDetails Ecliptic2Topocentric(double Lambda, double Beta, double Semidiameter, double Distance, double Epsilon, double Longitude, double Latitude, double Height, double JD)
        {
            double S = Globe.RhoSinThetaPrime(Latitude, Height);
            double C = Globe.RhoCosThetaPrime(Latitude, Height);

            //Convert to radians
            Lambda = CoordinateTransformation.DegreesToRadians(Lambda);
            Beta = CoordinateTransformation.DegreesToRadians(Beta);
            Epsilon = CoordinateTransformation.DegreesToRadians(Epsilon);
            Longitude = CoordinateTransformation.DegreesToRadians(Longitude);
            Latitude = CoordinateTransformation.DegreesToRadians(Latitude);
            Semidiameter = CoordinateTransformation.DegreesToRadians(Semidiameter);
            double sine = Math.Sin(Epsilon);
            double cose = Math.Cos(Epsilon);
            double cosBeta = Math.Cos(Beta);
            double sinBeta = Math.Sin(Beta);

            //Calculate the Sidereal time
            double theta = Sidereal.ApparentGreenwichSiderealTime(JD);
            theta = CoordinateTransformation.HoursToRadians(theta);
            double sintheta = Math.Sin(theta);

            //Calculate the Parallax
            double pi = Math.Asin(g_AAParallax_C1 / Distance);
            double sinpi = Math.Sin(pi);

            double N = Math.Cos(Lambda) * cosBeta - C * sinpi * Math.Cos(theta);

            TopocentricEclipticDetails Topocentric = new TopocentricEclipticDetails();
            Topocentric.Lambda = Math.Atan2(Math.Sin(Lambda) * cosBeta - sinpi * (S * sine + C * cose * sintheta), N);
            double cosTopocentricLambda = Math.Cos(Topocentric.Lambda);
            Topocentric.Beta = Math.Atan(cosTopocentricLambda * (sinBeta - sinpi * (S * cose - C * sine * sintheta)) / N);
            Topocentric.Semidiameter = Math.Asin(cosTopocentricLambda * Math.Cos(Topocentric.Beta) * Math.Sin(Semidiameter) / N);

            //Convert back to degrees
            Topocentric.Semidiameter = CoordinateTransformation.RadiansToDegrees(Topocentric.Semidiameter);
            Topocentric.Lambda = CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(Topocentric.Lambda));
            Topocentric.Beta = CoordinateTransformation.RadiansToDegrees(Topocentric.Beta);

            return Topocentric;
        }

    }
}
